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I. INTRODUCTION 

Given an expression for the total energy of an atom, 
molecule, or solid, 

Etot = T s + E e n + Eft + Fxc + (1) 

where the terms on the right-hand side represent the non¬ 
interacting kinetic, electron-nucleus, Hartree, exchange- 
correlation, and nucleus-nucleus energies, respectively, 
the search for the Slater determinant which minimizes 
E tot leads to one-electron Schrodinger equations 

(Mv 2 + v en (r) + v K (r) + t) xc (r)^ V’iO) = £»^»(r) (2) 

for the orbitals ifi and their energies £$. [For ease of no¬ 
tation, all formulas are given in spin-unpolarized form 
and for non-zero gap systems. N will denote the num¬ 
ber of (doubly) occupied orbitals.] In the Kohn-Sham 
(KS) 1 version of density functional theory (DFT)Pthe 
exchange-correlation potential v xc is calculated as the 
functional derivative of E xc with respect to the elec¬ 
tron density p (v xc = SE XC /Sp) and, as a consequence, 
v xc is a multiplicative potential (v xc ifi = v xc ifi), i.e., 
it is the same for all orbitals. Instead, in the gen¬ 
eralized KS (gKS) framework, formally introduced in 
Ref. 13, the derivative of E xc is taken with respect to 
ifi (?)xc ipi = $E XC as in the Hartree-Fock (HF) 
method. In other words, in the KS method the min¬ 
imization of the total energy [Eq. Q] is done with the 
constraint that the orbitals forming the (single) Slater de¬ 
terminant are solutions to a Schrodinger equation with a 
multiplicative potential, whereas in the gKS scheme this 
constraint is dropped. For exchange-correlation function¬ 
als E xc which depend explicitly only on p (and even¬ 
tually its derivatives) like in the local density approx¬ 
imation (LDA) 1 or generalized gradient approximation 


(GGA)P® the gKS method leads to the same multiplica¬ 
tive potential v xc as the KS method. However, for an 
orbital-dependent functional F xc , i.e., a functional which 
depends on the orbitals ^ not only via p, such as meta- 
GGA (see Ref. Eland references therein), self-interaction 
corrected^ or hybrid 8 functionals, the gKS method leads 
to a non-multiplicative (i.e., orbital-dependent) potential 
v xc — t’xc ,i as in the HF method. For such functionals, 
the calculation of SE XC /Sp , as required in the KS for¬ 
malism, is highly non-trivial, but possible by solving the 
optimized effective potential (OEP) equation. 9 

The focus of the present work will be on the multiplica¬ 
tive exchange potential v x . More specifically, approxi¬ 
mate semilocal exchange potentials will be compared to 
the exact exchange (EXX) potential obtained by means 
of the OEP method (called EXX-OEP thereafter), which 
has been implemented very recently 10 within the lin¬ 
earized augmented plane-wave 1 ® 5 ^ (LAPW) method for 
solids. 

The advantage of semilocal potentials, which depend 
on the local quantities p, Vp, V 2 p, or the kinetic-energy 
density t = J2iLi •V'lpi is that they are rather sim¬ 
ple to implement and lead to calculations which are much 
faster than EXX-OEP calculations or HF/hybrid calcu¬ 
lations with a non-multiplicative potential. The LDA 
exchange potential is for example a simple function of p, 
whereas in the case of GGA functionals, the correspond¬ 
ing potential becomes a function of p and its first two 
derivatives. 

The problem with the LDA and vast majority of GGA 
exchange functionals E x is that their functional deriva¬ 
tive v x barely resembles the EXX-OEP potential.^ 17 
Therefore several studies have focused on the search for 
better semilocal approximations for v x rather than E x . 
Among these studies there are the early works of Engel 
and Voskop 7 Baerends and co-workers 21 and the more 
recent works of Becke and Johnson,^ Staroverov and co- 
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workers pMSD Armiento et al. j 27 l 28 l and othersP^ES] 

It is worth mentioning that all these approximations 
for v x can be categorized in one of these two groups, 
namely, those which are functional derivative of a func¬ 
tional E x and those which are not (such potentials were 
termed stray in Ref. EU). Examples of exchange poten¬ 
tials which were modelled with the constrained to be a 
functional derivative are the ones from Engel and Vosko 17 
(EV93) and Armiento and KiimmeP^ (AK13), while the 
potentials from van Leeuwen and BaerendJ^ (LB94) and 
Becke and Johnson 22 (BJ) are stray potentials. Not con¬ 
straining a potential to be a functional derivative means 
much more freedom for its analytical form, however it has 
been shown that stray potentials have undesirable fea¬ 
tures both at the fundamental and practical level. 31 33 
Furthermore, attempts to turn a stray potential into a 
functional derivative without loosing too much of its orig¬ 
inal features have been rather unsuccessful up to now (see 
Refs. HQ E3 and l33|) . 

As already mentioned above, various semilocal ex¬ 
change potentials u x will be studied and compared to 
the EXX-OEP which will serve as reference. We will fo¬ 
cus in particular on the BJ potential, which has been 
shown to reproduce quite well the EXX-OEP poten tial 
in atoms 22 23 and has been applied to molecule d 34 * 35 * and 
solidsa s well as modified to improve the results in var¬ 
ious casesl22I2ZLiZI 4 l] Furthermore, in an attempt to be as 
close as possible to the EXX-OEP potential, a more gen¬ 
eral form of the BJ potential will be proposed. 

The paper is organized as follows. Section [IT] gives a 
summary of the EXX-OEP method and introduces the 
tested semilocal exchange potentials, while the computa¬ 
tional details are given in Sec. |III| Then, the results are 
presented and discussed in Sec. |IV| Finally, Sec. [V| gives 
the summary and an outlook for possible improvements. 


II. THEORY 


A. Optimized effective potential 


As mentioned in the Introduction, the calculation of 
the multiplicative exchange-correlation potential u xc = 
SE^ C /Sp for any orbital-dependent functional E^ c can 
be achieved by solving the integro-differential OEP 
equation! 9 for u xc (see Ref.|42]for a review), which is given 
in general by 


/x(r,r>xc(r')dV = A xc (r), 


with 


A, 


=<') = ?[/ 


<5£ xc fo/>j(r') 
r') <5veff(r) 

SE XC $Ei 


c.c. a r 


5si <5v e ff(r)J ’ 


( 3 ) 


d 3 -' 


( 4 ) 


where v c {[ = v en + vh + 't’xc is the KS effective potential 
and 


X(r,r') 


fy(r) 

^efr(r') 

0 A ^ (r)ipj (r)ipj (r')ipi(r') 

2 X X -777-+c-c<5) 


is the KS (non-interacting) density response function. 

So far, the OEP method has been applied mostly to 
the EXX energy 


N N 


j? xx = -E£// 


(r')V’i(r') , 3 , 3 / 

-:- — - a rd r , 


=i j=i ■ 


(6) 

which has the same analytic form as the HF exchange 
energy, but is evaluated with the KS orbitals instead of 
the HF orbitals. In the case of EXX, the sum over i in 
Eq. @ runs over the occupied orbitals only [Eq. § does 
not depend on unoccupied orbitals] and SE^ C /Ssi = 0. 


Talman and ShadwiclP^ were the first who reported 
EXX-OEP calculations (on spherical atoms). Initially, 
EXX-OEP was proposed as an approximation to HF 
in order to get rid of the non-multiplicative HF poten¬ 
tial. Later, it was recognized that the EXX-OEP method 
represents also the exact exchange method within the 
KS DFT framework (see, e.g., Ref. 16 and references 
therein). Since then EXX-OEP has attracted more and 
more attention. For solids the first EXX-OEP imple¬ 
mentation was reported by Kotani.® Subsequent reports 
of OEP calculations on solids (the focus of the present 
work) can be found in Refs. IT0TH21145H6T1 


The implementation of the OEP equation is rather 
complicated and its solution prone to instabilities in par¬ 
ticular if localized basis functions are used. 62 64 Recently, 
the implementation of the EXX-OEP method within the 
LAPW method has been reported.^ 10 12 ^ It employs an 
auxiliary basis, the mixed product basis, for representing 
the OEP. As discussed in detail in Refs. fT0landl65l in or¬ 
der to obtain a stable and physical EXX-OEP potential, 
the orbital (LAPW) and auxiliary (mixed product) basis 
sets have to satisfy a basis-set balance condition. This 
condition is fulfilled when the orbital basis set is con¬ 
verged with respect to the auxilary basis set and usually 
demands large orbital basis sets. The usage of (uneco- 
nomically) large LAPW basis sets can be avoided if the 
response of the LAPW basis functions is explicitly taken 
into account in the calculation of the KS orbital response 
S'lpi / in Eq. (El) and the calculation of the KS density 
response [Eq. ([5f]. 41 12 In this way, a much faster con¬ 
vergence is achieved with respect to basis set size (and 
number of unoccupied states) and the basis balance con¬ 
dition is fulfilled with smaller orbital basis sets. 
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B. Semilocal potentials 

Among the considered semilocal exchange potentials, 
LDAS as well as the GGAs B88® PEEP EV93P and 
AKlPl are functional derivatives of exchange-energy 
functionals that have the generic form 

Ex = (f) / p4/3 W Fx ( s ( r )) rf3 G ( 7 ) 

where F x (s) is the so-called exchange enhancement fac¬ 
tor which depends on the reduced density gradient s = 

|Vp| / ^2 (37r 2 ) 1 / 3 /? 4 / 3 ^. LDA is the exact form for the 

homogeneous electron gas and corresponds to F x (s) = 1. 
As shown in Fig. [TJ the enhancement factors of the GGA 
functionals are larger than one, thus correcting the ten¬ 
dency of LDA to underestimate the magnitude of the 
exchange energy. Compared to the standard PBE func¬ 
tional, EV93 and AK13 are much stronger. Note that at 
s = 0 all factors F x (s) reduce to one in order to satisfy 
the homogeneous electron gas limit given by LDA. B88 
and PBE, which are among the most popular GGA func¬ 
tionals for calculating the properties of molecules and 
solids, respectively, were constructed without consider¬ 
ing the quality of the potential v x . For EV93 and AK13, 
however, the emphasis was put on v x . The parameters in 
EV93 were determined by a fit to EXX-OEP potentials in 
atoms pl while Armiento and Kiimmel were able to find 
an analytical form for AK13 such that u x changes discon- 
tinuously at integer particle numbers.^ Both EV93 and 
AK13 were shown to improve over the standard LDA and 
PBE functionals for the band gaps in sohdsP^ ^ 6 7 

In addition to these potentials, we consider in this work 
the BJ potential 22 which is of stray type (see Refs. 02] and 
l37j) and has the form 



FIG. 1. (Color online) The enhancement factors F x (s) [see 
Eq. Q] of the different exchange functionals considered in 
this work. 


with 


D( r) = t( r) - 


1 |VKr)| 2 

8 p( r) 


After x is calculated, b in Eq. ( fTo] ) is given by 


b( r) 


x 3 (r)e 

47rp(r) 


1/3 


(13) 
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or the Becke-Roussel potential^ 

0 -x(r) _ i 






x(r)e 


-x(r) 


(9) 


( 10 ) 


The function x in Eq. (10) is calculated by solving (at 
each point of space) the nonlinear equation (or using the 
analytical interpolation formula for x from Ref. HO]) 


c(r)e 2a 4 r )/ 3 1 /tj- \ 2/3 p 5 / 3 ( r ) 

_ 3 V2 ) 


x 


( r )-2 


Q(r) 


where 


Q( r ) = 3 ^ (V 2 p(r) - 4rfD(r)) 


( 11 ) 

( 12 ) 


The parameter 7 in Eq. © has to be set to 1 or 
0.8 in order to recover the exact exchange potential of 
the hydrogen atom or the homogeneous electron gas, 
respectively.^ Note that since the BR potential and the 
second term of Eq. depend on the kinetic-energy den¬ 
sity £, they are of the semilocal meta-GGA form, 7 ^ while 
the Slater potential [Eq. (|9|] is nonlocal in the sense that 
the calculation of v 8 at r requires the value of quanti¬ 
ties (the occupied orbitals) at all points of space r'. For 
closed-shell atoms it was shown t hat the BR potential is 
very close to the Slater potential . 22 69 In the rest of this 
work, we focus on the BJ-based potentials using the BR 
potential. 

Severa l mod ifications of the BJ potential have been 
proposed! 23 * 27 * 38 * 3 -^ For instance, in Ref. 38] we proposed 
the modified BJ (mBJ) potential 


< BJ (r) = c BR (r) + (3c - 2) , (15) 

where c is a parameter that was introduced to improve 
the agreement with experiment for the band gaps in 
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solids and that was parameterized using the average of 
|Vp| / p in the unit cell. 

As pointed out by Rasanen et al. in Ref. [39] the BJ 
potential is not gauge-invariant, does not show the cor¬ 
rect asymptotic behavior at r —» oo in finite systems, and 
the correction to the Slater (or BR) term is not zero for 
one-electron systems as it should be. In order to cure 
these deficiencies, they proposed an universal correction 
(UC) to the BJ potential. For systems with zero current 
density J (like those considered in this work), the UC 
consists of replacing t by D [Eq. ( p~3] )] in the second term 
of Eq. ([ 8 |. 

In an attempt to propose in the present work a semilo¬ 
cal exchange potential which can reproduce accurately 
the EXX-OEP, we will consider a more general form of 
the BJ and mBJ potentials, called generalized BJ (gBJ) 
thereafter: 


J ( = 


(r) = ci£ K (r)+(3c-2)- 


1(f) 


3\l/3 


t p (r) 


(&(3tt 2 ) 2/3 )V 5 - 1)/3 (r)’ 

(16) 

which, in addition to the two parameters 7 [in 
Eq. (12)] and c as in mBJ, contains a third one (p) whose 
value is 0.5 in BJ and mBJ. The form of the second term 
of Eq. (16) was chosen such that the LDA exchange po¬ 
tential is recovered for constant electron densities (and 
if 7 = 0.8, see above). As a modification of Eq. (16), 
we will also consider its variant where t in the second 
term is replaced by D (UC, Ref. 39), leading to the gB- 
JUC potential. The parameters 7 , c, and p of the gBJ 
and gBJUC potentials were varied around the standard 
BJ values within the intervals [0.4,1.4], [1.0,1.4], and 
[0.35, 0.65] and in steps of 0.2, 0.1, and 0.05, respectively. 
In the following, BJ will denote the unmodified potential 
given by Eq. (| 8 | using BR with 7 = 0.8 and similarly 
for BJUC. The values of the parameters in the gBJ and 
gBJUC potentials will be specified in this order: ( 7 , c,p). 

Regarding the influence of the parameters 7 , c, and 
p on the shape of the gBJ potential, we have generally 
observed that an increase of one or another of the pa¬ 
rameters leads to more pronounced variations, and this 
effect is rather similar for the three parameters (see Fig. [ 2 ] 
for an illustrative example in the diamond phase of Cp 
Nevertheless, by looking more closely at Fig. [ 2 | we can 
notice some differences in the way the parameters 7 , c, 
and p modify the potential. For instance, when 7 is in¬ 
creased [Fig. [^a)], the intershell peak at d ~ 0.4 A gets 
more spiky, while the value of c affects the potential in 
a broader region of space [Fig. [2|b)]. This is the reason 
why we have found it useful to use three parameters in¬ 
stead of only one or two in our attempt for reproducing at 
best the EXX-OEP results with the gBJ(UC) potential. 

The effect of the UC is shown in Fig. |3][a) by com¬ 
paring the BJ and BJUC potentials in C. Rather large 
differences between the two potentials are visible in the 
core region (d < 0.4 A in this example) where the BJUC 
potential is much more attractive than BJ. As a conse¬ 
quence, the core states are bound stronger when the UC 



FIG. 2. (Color online) gBJ( 7 , c,p) exchange potentials in C 
plotted starting at a distance of 0.2 A from the atom at site 
(1/4,1/4,1/4) in the [111] direction. The center of the unit 
cell is at d = 2.32 A. The potentials were shifted such that 
fee u Mr)d 3 r = 0. 


is used. We made the same observation for all other in¬ 
vestigated solids (see also Fig. 3 of Ref. [39] for the Ne 
atom). In Fig. [3](b) , the kinetic-energy density t and D 
[Eq. p| )] are compared by showing the ratio t/D, where 
we can see that it is indeed in the core region that t/D 
differs the most from 1 . Actually, the term |Vp | 2 / ( 8 p) in 
D is the von Weizsacker 72 kinetic-energy density which 
is equal to the exact kinetic-energy density t in regions 
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FIG. 3. (Color online) (a) BJ and BJUC exchange potentials 
in C. (b) Ratio t/D for C. The path and potential shift are 
as in Fig. [l] 

of space dominated by a single orbital. Therefore, in 
such regions D = 0 and the BJUC potential reduces to 
the BR (or Slater) potential which is by far too negative 
compared to EXX-OEP.^ Note that the (smaller) differ¬ 
ences between t and D in the region 0.4—1.4 A also affect 
the potentials, however, due to the alignment of the dif¬ 
ferent potentials (f n v x (r)d 3 r = 0), the differences be¬ 
tween BJ and BJUC are transferred to d > 1.2 A. 


III. COMPUTATIONAL DETAILS 

The calculations with the semiloc al p otentials and 
EXX-OEP were done with the WIEN2fPand FLEURP 
codes, respectively. Since the two codes use the same 
basis set (LAPW^^), it was also possible to calculate 
the EXX-OEP orbitals with WIEN2k by fixing the po¬ 
tential u x to the EXX-OEP read from a file (contain¬ 
ing the radial functions and Fourier coefficients of the 
spherical harmonics and plane-wave expansions of the 
potential) generated by FLEUR. Despite some (small) 
technical differences between the two codes and different 
computational parameters used for the calculations (e.g., 
basis sets), we observed for all solids, that running a PBE 


calculation with WIEN2k as usual or with the FLEUR- 
generated PBE potential leads to very close results (e.g., 
the transition energies differ by less than 0.02 eV). There¬ 
fore, we are convinced that this procedure of calculating 
orbitals using a potential generated from another code 
is reliable in terms of accuracy. In addition, HF calcu¬ 
lations with the WIEN2k code were also done. Note, 
however, that in the current implementation of the HF 
method) 75 -* the core electrons experience a semilocal po¬ 
tential, like in other impl ement ations of the HF method 
with the LAPW basis setP 6 77 The k-mesh for the inte¬ 
gration of the Brillouin zone and size of the basis sets 
were chosen to be converged for the purpose of our work. 

The set of solids that we will consider consists of the 
nonmagnetic cubic (the space group, number of atoms 
in the primitive unit cell, and cubic lattice constant are 
indicated in parenthesis) C (Fd3m, two atoms, 3.57 A), 
Si (Fd3m, two atoms, 5.43 A), MgO (Fm3m, two atoms, 
4.23 A), BN (F43m, two atoms, 3.62 A), and CU 2 O 
(Pn3m, six atoms, 4.27 A). C, Si, and BN are cova¬ 
lent, while MgO and CU 2 O are ionic. Also included in 
our test set is NiO whose type-II antiferromagnetic order 
along the [111] direction reduces the symmetry from cu¬ 
bic (Fm3m, two atoms, 4.17 A) to rhombohedral (R3m, 
four atoms). All these solids are nonmetallic and while 
most of them are simple sp-type semiconductors or in¬ 
sulators, two of them, namely CU 2 O and NiO, represent 
more stringent tests. 

NiO is a rather difficult system to describe 
theoretically 78 since the Ni-3d electrons are strongly 
correlated as it is generally the case in magnetic 3d- 
transition-metal oxides. Due to their inherent self¬ 
interaction errorf 7 the semilocal functionals perform par¬ 
ticularly bad in such systems and more advanced meth¬ 
ods like DFT+U 7 ^are commonly used. In Refs. [75] and 
1501 we showed that a correct description of the band gap 
and electric-field gradient (EFG) in CU 2 O could only be 
achieved with the hybrid functionals, while the results 
obtained with the LDA, GGA, LDA-j-U, and mBJ meth¬ 
ods were qualitatively wrong. 


IV. RESULTS AND DISCUSSION 

We quantify the difference between the semilocal ex¬ 
change potentials and the reference EXX-OEP by com¬ 
paring the EXX total energy and the electronic structure 
for our test set of solids. The electronic structure of the 
solids is assessed in terms of the band transition across 
the band gap, the position of the core electrons, and the 
density of states (DOS). Furthermore, as a measure of 
the similarity of the electron density we compare the re¬ 
sulting EFG in CU 2 O and the magnetic moment in NiO. 
We start with the discussion of the EXX total energy. 
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TABLE I. EXX total energy E t ^ xx (in Ry/cell) calculated with orbitals generated from various exchange potentials. The values 
for the semilocal potentials are the differences with respect to EXX-OEP, and a positive value indicates that EXX-OEP leads 
to a more negative energy as it always should. 


Potential 



c 

Si 

BN 

MgO 

Cu 2 0 

NiO 

EXX-OEP 



-151.592 

-1158.353 

-158.623 

-550.129 

-13527.744 

-6377.723 

LDA 



0.042 

0.079 

0.047 

0.079 

0.576 

0.949 

PBE 



0.026 

0.040 

0.027 

0.037 

0.351 

0.619 

B88 



0.026 

0.040 

0.026 

0.036 

0.357 

0.620 

EV93 



0.017 

0.015 

0.015 

0.010 

0.206 

0.420 

AK13 



0.029 

0.027 

0.023 

0.014 

0.152 

0.312 

BJ 



0.008 

0.019 

0.008 

0.015 

0.177 

0.395 

BJUC 



0.074 

0.096 

0.067 

0.072 

0.256 

0.642 

gBJ(0.6,1.0,0.60; 

a 


0.003 

0.000 

0.002 

0.001 

0.154 

0.264 

gBJ(1.4,1.1,0.50; 

b 


0.014 

0.054 

0.013 

0.037 

0.286 

0.257 

gBJ(0.4,1.3,0.65; 

c 


0.159 

0.240 

0.148 

0.199 

0.726 

0.335 

gBJUC(1.4, 1.2,0.50^J 

0.202 

0.272 

0.202 

0.257 

0.786 

0.757 


a Good compromise for the EXX total energy of C, Si, BN, MgO, and CU 2 O. 
b Good compromise for transition energies in C, Si, BN, and MgO. 
c Good compromise for transition energies and Ni magnetic moment in NiO. 
d Good compromise for transition energies and EFG in CU 2 O. 

A. EXX total energy 

The EXX total energy E^ xx is calculated with orbitals 
either generated from the multiplicative EXX-OEP or 
the semilocal exchange potentials. The obtained total 
energies are shown in Table |TJ The lowest EXX total en¬ 
ergy is obtained by using the EXX-OEP orbitals, which 
was expected since the EXX-OEP is also the solution of 
the equation 5Ef xx / = 0.® The LDA orbitals lead 

to energies which are higher by 0.04-0.08 Ry for C, Si, 

BN, and MgO, 0.6 Ry for CU 2 O, and 0.9 Ry for NiO. 

All sets of GGA (PBE, B 88 , EV93, and AK13) orbitals 
improve upon LDA by reducing the difference with re¬ 
spect to EXX-OEP by a factor of two up to four. On 
average EV93 and AK13 yield total energies which are 
closer to the EXX-OEP energy than PBE and B 88 . The 
BJ potential [Eq. (J 8 |] shows a rather similar performance 
as EV93 and AK13, while BJUC leads to total energies 
that are sometimes even worse than LDA. 

The results for the gBJ potential [Eq. ( p~6] )] are shown 
for a few selected sets of parameters ( 7 T~p). With the 
parameters (0.6,1.0,0.60) the results are close to optimal 
(within the space of parameters) for E^ xx and all solids 
except NiO for which an increase of c to 1.2 or 1.3 would 
further reduce the difference with respect to EXX-OEP 
by a factor of two. It should be stressed that the er¬ 
ror obtained with gBJ(0.6,1.0,0.60) is only of the order 
of 0.001%, i.e., below 1-3 mRy for the light solids with¬ 
out transition-metal atoms. Nevertheless, as shown be¬ 
low, such a good agreement for the total energy does not 
necessarily mean a good agreement with EXX-OEP for 
other quantities like the transition energies, which require 
other sets of parameters (q,c,p) (also shown in Table [I]). 


It should be also mentioned that showing the results for 
the parameters ( 0 . 6 , 1 . 0 , 0.60) is only one choice among a 
few others which lead to a similar (albeit maybe slightly 
worse overall) agreement with EXX-OEP. For instance, 
by varying only c with respect to the original BJ po¬ 
tential (i.e., considering mBJ) the results for E^ xx are 
also very good with c = 1.1. The gBJUC orbitals lead 
systematically to very high EXX total energies whatever 
the parameters (q,c,p) are. Actually, this is related to 
the poor reproduction of the EXX-OEP potential by gB¬ 
JUC in the region close to the nuclei (see below) which 
substantially affects the total energy. 


B. Electronic structure 

We now turn to the discussion of the electronic struc¬ 
ture and focus first on the comparison of the direct tran¬ 
sition energies across the band gap. 


1 . Transition energies 

For each solid the direct transition energies are calcu¬ 
lated at three k-points in the Brillouin zone (expressed 
in primitive basis for NiO and conventional basis for 
the other solids): T = (0,0,0), X = (0,1,0), and 
L = ( 1 / 2 , 1 / 2 , 1 / 2 ) for C, Si, BN, and MgO. T = (0, 0, 0), 
X = ( 0 , 1 / 2 , 0 ), and M = ( 1 / 2 , 1 / 2 , 0 ) for Cu 2 0. T = 
( 0 , 0 , 0 ), L = ( 0 , 1 / 2 , 0 ), and F = ( 0 , 1 / 2 , 1 / 2 ) for NiO. 
The mean error (ME) and mean absolute error (MAE) 
with respect to the EXX-OEP is shown in Table [II] for 
the different solids and potentials. Applying the LDA the 
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TABLE II. Statistics on direct transition energies A^k = ejy+i — £jv at three different k-points (specified in the text). The 
values for EXX-OEP are the mean over the three k-points of the transition energy (J^ k Aejp 00 ' OEP ), while for the semilocal 
potentials the values are the MAhj^jand Mhj^Jwith respect to EXX-OEP. All values are in eV. 


Potential 



C 


Si 


BN 

MgO 

Cu 2 0 

NiO 


EXX-OEP 


9.80 

3.49 

10.93 

9.50 

3.08 

5.37 


MAE 

ME 

MAE 

ME 

MAE 

ME 

MAE 

ME 

MAE 

ME 

MAE 

ME 

LDA 


0.62 

-0.62 

0.68 

-0.68 

0.97 

-0.97 

1.95 

-1.95 

1.23 

-1.23 

3.68 - 

-3.68 

PBE 


0.32 

-0.32 

0.41 

-0.41 

0.62 

-0.62 

1.30 

-1.30 

1.05 

-1.05 

3.11 

-3.11 

B88 


0.32 

-0.32 

0.40 

-0.40 

0.60 

-0.60 

1.26 

-1.26 

1.04 

-1.04 

3.08 

-3.08 

EV93 


0.31 

-0.18 

0.23 

-0.17 

0.41 

-0.41 

0.81 

-0.81 

0.98 

-0.98 

2.78 

-2.78 

AK13 


0.30 

-0.03 

0.38 

0.18 

0.27 

-0.11 

0.19 

0.19 

0.82 

-0.82 

2.27 

-2.27 

BJ 


0.27 

-0.27 

0.38 

-0.38 

0.45 

-0.45 

1.01 

-1.01 

0.95 

-0.95 

2.53 

-2.53 

BJUC 


0.38 

-0.38 

0.53 

-0.53 

0.45 

-0.45 

0.44 

-0.44 

0.42 

-0.42 

2.64 

-2.64 

gBJ(0.6,1.0,0.60; 

C 

0.13 

-0.13 

0.17 

-0.17 

0.29 

-0.29 

0.71 

-0.71 

1.01 

-1.01 

2.22 

-2.22 

gBJ(1.4,1.1,0.50; 

3 

0.14 

0.14 

0.06 

0.02 

0.07 

0.07 

0.19 

-0.19 

0.87 

-0.87 

1.49 

-1.35 

gBJ(0.4,1.3,0.65; 

j 

0.86 

0.86 

1.42 

1.42 

1.13 

1.13 

1.67 

1.67 

1.19 

-1.19 

0.28 

0.11 

gBJUC(1.4,1.2, 0.50^j 

0.10 

0.10 

0.01 

-0.01 

0.30 

0.30 

0.88 

0.88 

0.06 

0.06 

1.70 

-1.54 


a MAE = |A£^ emilocal - AsE xx -° ep |. 
b ME = (A£^ emilocal - A£ EXX_OEP ). 

c Good compromise for the EXX total energy of C, Si, BN, MgO, and CU 2 O. 
d Good compromise for transition energies in C, Si, BN, and MgO. 
e Good compromise for transition energies and Ni magnetic moment in NiO. 
f Good compromise for transition energies and EFG in CU 2 O. 


MAE is in the range of 0.6-3.7 eV where the largest error 
is for NiO. Actually, it is well known 49 that LDA strongly 
underestimates the band gap with respect to experiment 
and EXX-OEP. The GGA, and in particular EV93 and 
AK13, improve over LDA by reducing the MAE by a 
few 0.1 eV for C, Si, BN, and CU 2 O or more than 1 eV 
for MgO and NiO, but overall the errors remain rather 
substantial. BJ and BJUC perform similarly to EV93 or 
AK13. For all these potentials except AK13, the nega¬ 
tive sign of the ME and its magnitude which is equal to 
the MAE in most cases indicate that the deviation from 
EXX-OEP corresponds to an systematic underestimation 
of the transition energies. 

For gBJ, we found that the parameters (1.4,1.1,0.50) 
lead to a very small MAE (below 0.2 eV) for C, Si, BN, 
and MgO. For CU 2 O and NiO different sets of parame¬ 
ters are required. For CU 2 O it was not possible to find 
a combination of the parameters (within the considered 
ranges) that reduces the MAE below 0.7 eV. However, 
by considering the gBJ potential with the UC (gBJUC), 
we were able to improve substantially the results for the 
transitions energies. For instance (see Table [TI|, with 
the parameters (1.4,1.2,0.50), gBJUC leads to a MAE 
of 0.06 eV for the transition energies of CU 2 O. 

In the case of NiO, a substantial improvement for the 
transition energies can be obtained if the parameter c is 
increased to at least 1.2. For instance, with the param¬ 
eters (0.4,1.3,0.65) the MAE on the transition energies 
is below 0.3 eV, which is one order of magnitude smaller 
than with the other methods. We mention that for NiO, 
the values of the parameters 7 and p have little influence 


on the results and only an increase of c can lead to a clear 
improvement. 


2. Core states 

We proceed by discussing the effect of the different 
potentials on the binding energies of the core electrons. 
Table m shows the averaged energetic position of the 
core states with respect to the Fermi energy for the dif¬ 
ferent solids and potentials. For the definition of the 
mean absolute relative error (MAR E) a nd mean relative 
error (MRE) see caption of Table |III[ As observed for 
the EXX total energy and the transition energies, all 
GGA exchange potentials improve over LDA by reduc¬ 
ing the MARE below 0.5% for most solids. The posi¬ 
tive MRE for LDA, PBE, and B 88 indicate that these 
potentials lead to core states which are typically bound 
to loosely with respect to EXX-OEP. For EV93, AK13, 
and BJ there is no systematic trend. Among the four 
selected parameterizations of the gBJ potential it turns 
out that (0.4,1.3,0.65) (optimized for NiO) leads overall 
to a rather clear improvement over the LDA and GGA 
potentials. The accuracy obtained with the set of param¬ 
eters (0.6,1.0, 0.60) (optimized for the EXX total energy) 
is satisfying except for C. The results obtained with the 
gBJUC-based potentials are extremely inaccurate, which 
is, as already mentioned above, due to the very poor re¬ 
production of the EXX-OEP close to the nuclei, leading 
to core states that are too low in energy. 















TABLE III. MARE and MREp](with respect to EXX-OEP and in %) for the energy position of the core states with respect 
to the valence band maximum (A^core,i = Score,i — £vbm)- The MARE and MRE are over all core states in the solid: C (Is), 
Si (Is), BN (B: Is; N: Is), MgO (Mg: Is; O: Is), Cu 2 0 (Cu: Is, 2s, 2p; O: Is), NiO (Ni: Is, 2s, 2p; O: Is). A negative MRE 
means that on average the core states are deeper in energy than EXX-OEP. 


Potential 


c 


Si 


BN 


MgO 

Cu 2 0 

NiO 

MARE 

MRE 

MARE 

MRE 

MARE 

MRE 

MARE 

MRE 

MARE 

MRE 

MARE 

MRE 

LDA 


1.67 

1.67 

0.83 

0.83 

2.46 

2.46 

1.28 

1.28 

0.70 

0.70 

0.72 

0.70 

PBE 


0.28 

0.28 

0.31 

0.31 

0.95 

0.95 

0.49 

0.49 

0.37 

0.37 

0.52 

0.45 

B88 


0.20 

0.20 

0.28 

0.28 

0.85 

0.85 

0.45 

0.45 

0.34 

0.34 

0.51 

0.42 

EV93 


0.27 

-0.27 

0.21 

0.21 

0.54 

0.40 

0.27 

0.26 

0.35 

0.32 

0.48 

0.43 

AK13 


1.25 

-1.25 

0.07 

-0.07 

0.66 

-0.66 

0.36 

-0.19 

0.39 

0.12 

0.50 

0.28 

BJ 


1.13 

-1.13 

0.11 

-0.11 

0.45 

-0.42 

0.46 

-0.24 

0.27 

-0.04 

0.53 

0.08 

BJUC 


7.82 

-7.82 

2.08 

-2.08 

7.44 

-7.44 

3.56 

-3.56 

1.50 

-1.50 

1.17 

-1.17 

gBJ(0.6,1.0,0.60; 


0.74 

-0.74 

0.03 

-0.03 

0.56 

0.00 

0.39 

-0.05 

0.20 

-0.02 

0.48 

0.08 

gBJ(1.4,1.1,0.50; 


1.61 

-1.61 

0.06 

-0.06 

1.03 

-1.03 

0.56 

-0.40 

0.64 

0.20 

0.63 

0.43 

gBJ(0.4,1.3,0.65; 

eT 

0.11 

-0.11 

0.04 

0.04 

0.67 

0.62 

0.24 

0.19 

0.17 

0.00 

0.40 

0.15 

gBJUC(1.4,1.2, 0.50^j 

12.96 

-12.96 

3.35 

-3.35 

12.96 - 

-12.96 

6.06 

-6.06 

2.21 

-2.21 

1.72 

-1.62 


a MARE = E-° re 100 cal - A^xx-oep| / | A£ exx-oep|. 

b MRE = E'° re 100 (^e S co^ cal - Ae=*£ OEP ) / | A e EX e X :° EP J- 
c Good compromise for the EXX total energy of C, Si, BN, MgO, and Cu 2 0. 
d Good compromise for transition energies in C, Si, BN, and MgO. 
e Good compromise for transition energies and Ni magnetic moment in NiO. 
f Good compromise for transition energies and EFG in Cu 2 0. 


3. Density of states 

The comparison of the electronic structure obtained 
with the different exchange potentials focussed so far on 
the transition energies and the core states. In order to as¬ 
sess the differences in the electronic structure on a wider 
energy spectrum of the valence states, we compare the 
density of states of NiO and CU 2 O around the Fermi en¬ 
ergy. We picked out these two solids from our test set, 
since the largest changes in the DOS with respect to the 
applied potential can be observed for these two solids. 
For C, Si, BN, and MgO the basic structure of the DOS 
remains very similar independent from the applied poten¬ 
tial (of course, apart from a rigid shift of the conduction 
bands). 

Figures [4] and [5] show the DOS of CU 2 O and NiO, re¬ 
spectively, for a few selected potentials. In the case of 
CU 2 O, we can clearly see that the gBJUC(1.4,1.2,0.50) 
potential (very good for the transition energies and EFG, 
see below) leads to the best agreement with EXX-OEP, 
which is particularly true for the partial Cu-3d DOS in 
the range from —3 to 0 eV below the Fermi energy. The 
Cu-3d DOS obtained with the other semilocal potentials, 
including gBJ without UC, are too broad by about 1 eV. 
gBJUC(1.4,1.2,0.50) also leads to correct positions of 
both the 0-2 p (extending from —7 to —5 eV) and the 
conduction band Cu-4s states. The HF DOS differs sig¬ 
nificantly from the other calculations employing a local 
exchange potential including EXX-OEP. It is well known 
that the HF method systematically leads to band gaps 


which are by far too large compared to experiment. In 
the case of CU 2 O, the HF band gap amounts to 10.4 eV, 
while it is only 2.17 eV in experiment 81 and 1.44 eV with 
EXX-OEP. In the occupied part of the spectrum, it is 
obvious that the main position of the Cu-3d peaks are 
much lower in energy (below —4 eV), while the bands 
in the energy range from —4 to 0 eV are a mixture of 
Cu-3d and 0-2p states. For other systems, a comparison 
between the HF and EXX-OEP occupied spectrum can 
be found in, e.g., Refs. 02j 82]. and 1551 

For NiO, the semilocal methods lead to DOS that differ 
markedly from the EXX-OEP DOS. As already discussed 
in Refs. H2j and 15011 he spin-up highest valence bands (be¬ 
tween —0.8 and 0 eV) and lowest conduction bands in the 
EXX-OEP DOS are of Ni-3d character coming from the 
Ni atom with more spin-down electrons (Ni2 with t\ g oc¬ 
cupied and ej emtpy), while the spin-up states between 
—7.5 and —2 eV are mixtures of Ni-3d from the Ni atom 
Nil (t\ g and ej fully occupied) and 0-2p. Therefore, 
EXX-OEP leads to a clear energy separation between 
the spin-up and spin-down Ni-3d states of the same Ni 
atom. In the PBE DOS the position of the conduction 
states is much too low and there is no Ni-3d peak similar 
to the one at ~ —7.5 eV in the EXX-OEP DOS. In ad¬ 
dition, there is very little energy separation between the 
spin-up 3d states coming from the two Ni atoms. The 
gBJ(0.4,1.3, 0.65) potential leads to a good band gap and 
a separation of ~ 0.5 eV between the spin-up 3d states 
from the two Ni atoms, however, there is still no Ni-3d 
peak at the lower part of the valence DOS, which can 
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Energy (eV) 

FIG. 4. (Color online) DOS of CU 2 O. The Fermi energy is set 
at zero. 


only be obtained by the LDA+l® or HF/hybricPHHH 
methods. We mention that the gBJ potentials with small 
values of c (1.0 or 1.1) and the gBJUC potentials do not 
produce any energy separation between the 3d states of 
the two Ni atoms. The valence part of the HF DOS starts 
at —10 eV and about five sharp Ni-3d peaks are equally 
distributed in the energy range —10 to —5 eV, while the 
DOS from —3 to 0 eV is exclusively of 0-2 p character. 
The HF band gap is 13.9 eV, which is in fair agreement 
with previous HF results!^*®* In experiment, a gap of 
4.0-4.3 eV is observed* 88 * 89 * while EXX-OEP gives rise to 



-10 -8 -6 -4 -2 0 2 4 6 8 


Energy (eV) 

FIG. 5. (Color online) Spin-up DOS of NiO. The Fermi energy 
is set at zero. 


3.54 eV. 


C. EFG of C112O and magnetic moment in NiO 


As shown in Refs. [94] and [55] the EFG is mainly de¬ 
termined by the non-spherical electron density close to 
the nucleus. Since the density of the core electrons is (by 
construction) purely spherical, it is the electron density 
of the valence states that determines the EFG. More¬ 
over, in the case of a 3d-transition metal like Cu, the 
valence electron density in a region of a few tenths of 
an Angstrom from the nucleus is decisive. Hence, by 
comparing the EFG of CU 2 O for the different poten¬ 
tials we indirectly measure the difference in the non- 
spherical part of the valence electron density. The results 
for the EFG of Cu are shown in Table [IV Similarly to 
the transition energies in CU 2 O (see Sec. IVB1), only 


the gBJ potential with the UC (gBJUC) is able to re- 
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TABLE IV. Spin magnetic moment pjg 1 (in fi b) inside the Ni 
atomic sphere of radius 1.016 A in NiO and EFG of Cu (in 
10 21 V/m 2 ) in Cu 2 0. 


Method 



// Ni 

MS 

Total 

EFGcu 

P-P 

d-d 

EXX-OEP 



1.91 

-17.7 

-25.0 

7.2 

LDA 



1.30 

-4.7 

-15.7 

10.8 

PBE 



1.43 

-5.6 

-16.3 

10.5 

B88 



1.43 

-5.6 

-16.3 

10.5 

EV93 



1.51 

-6.8 

-17.5 

10.4 

AK13 



1.58 

-8.1 

-18.5 

10.1 

BJ 



1.53 

-7.4 

-17.7 

10.2 

BJUC 



1.41 

-11.3 

-19.4 

7.9 

gBJ(0.6,1.0,0.60; 

a 


1.61 

-7.0 

-17.6 

10.5 

gBJ(1.4,1.1,0.50; 

b 


1.66 

-8.3 

-19.3 

10.8 

gBJ(0.4,1.3,0.65; 

c 


1.86 

-5.1 

-18.2 

13.0 

gBJUC(1.4,1.2,0.50 rl 

1.59 

-15.1 

-22.2 

6.8 

HF 



1.88 

-17.0 

-25.0 

7.9 

Expt. 



1.90 ±0-1] 

9.fl 




a Good compromise for the EXX total energy of C, Si, BN, MgO, 
and CU 2 O. 

b Good compromise for transition energies in C, Si, BN, and 
MgO. 

c Good compromise for transition energies and Ni magnetic 
moment in NiO. 

d Good compromise for transition energies and EFG in Cu 2 0. 

e Does not include the orbital contribution = 0.32 d= 0.05 fiB 
(Refs. l90l and I91D . 

f Only the magnitude is known. Calculated using the quadrupole 
moment Q ( 63 Cu) = 0.22 ( Refs. [921 and |93| . 


produce the EXX-OEP EFG qualitatively. For instance, 
with the parameters (1.4,1.2,0.50), gBJUC leads to an 
EFG of —15.1 x 10 21 V/m 2 , while for all other poten¬ 
tials (except BJUC), the magnitude of the EFG does 
not exceed 10 x 10 21 V/m 2 . For gBJUC(1.4,1.2, 0.50), 
not only the total EFG, but also its two main compo¬ 
nents (p-p and d-d) agree rather well with the EXX- 
OEP (and HF) values (see Table IV). A detailed anal¬ 
ysis of the UC will be provided in Sec. UVD] A calcula¬ 
tion with the non-multiplicative HF potential leads to an 
EFG of —17.0 x 10 21 V/m 2 which is relatively close to 
the EXX-OEP value and expected since the first-order 
change in the electron density due to the replacement 


..HF 


^EXX-OEP lg zer0 L 


The magnitude of the ex¬ 
perimental EFG amounts to 9.8 x 10 21 V/m 2 (Refs. 92 
and which is much smaller than the EXX-OEP or 
HF values. Thus, the impact of the electron correlation 
on the EFG is significant: the EXX-OEP value has to 
be reduced by the exact correlation functional nearly by 
its half. Furthermore, it was shown in Refs. cza Ea ebi- 
mm that the LDA, GGA, LDA +U, onsite-hyl)ri(ipIM2l 
and mBJ methods lead to an EFG in CU 2 O which is by 
far too small. The experimental value could only be ap¬ 
proached with full hybrid functionals. 



FIG. 6. (Color online) Two-dimensional plots of exchange 
potentials v x in a (110) plane of C. The potentials were shifted 
such that f cell v x (r)d 3 r = 0. The contour lines start at —2 Ry 
(blue color) and end at 1 Ry (red color) with an interval of 
0.2 Ry. 


Next we turn the discussion to the spin magnetic mo¬ 
ment /ig 1 of Ni in NiO (results in Table IV). In contrast to 
the EFG, fig 1 is determined by the difference of the spher¬ 
ical spin-up and -down electron densities in the atomic 
spheres. The best agreement with the EXX-OEP Ni spin 
magnetic moment is obtained by the gBJ potential with 
a c parameter of at least 1.2, which is in accordance with 
the observations for the EXX total energy and transition 
energies in Secs. |IV A| and |TVB 1 In fact, the parameters 
(0.4,1.3, 0.65) of the gBJ potential lead to fig 1 = 1.86 /ip, 
which is very close to the EXX-OEP, HF, and experimen¬ 
tal values. All other investigated potentials substantially 
underestimate the EXX-OEP spin-magnetic moment. 
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d (k) 


FIG. 7. (Color online) Exchange potentials v x in Si plotted 
from the vicinity of the atom at site (1/8,1/8,1/8) (d — 0 ) 
to (a) the center of the unit cell (d — 3.53 A) or (b) the mid¬ 
distance to the atom at site (5/8, 5/8,1/8) (d = 2.38 A). The 
potentials were shifted such that f cell v x (r)d 3 r = 0 . 


D. Analysis of the potentials 


In the following, the results of the previous subsections 
are set in relation to the spatial form of the different ex¬ 
change potentials. We start with discussing the exchange 
potential of C in the (110) plane (see Fig. [ 6 |. Since the 
LDA potential depends only on the electron density p, 
the corresponding contour plot is the most structure¬ 
less. In comparison to the other potentials it exhibits 
a more spherical shape around the C atoms and is less 
corrugated in the interstitial region. It is less attractive 
(i.e., negative) than EXX-OEP in the bonding region, 
but more attractive in the interstitial. Therefore, com¬ 
pared to LDA there is a transfer of electrons from the 
interstitial to the bonding region with EXX-OEP (see 
Sec. |TVEl ). The GGA potentials (PBE, B 88 , EV93, and 
AK13), that depend additionally on the first and second 
derivatives of p, show stronger spatial variations. For ex¬ 


ample, the PBE potential is more undulated than LDA, 
but features seen in the EXX-OEP are still reproduced 
too weakly or completely missing. (B 88 is not shown 
since its contour plot is very similar to the PBE plot.) 
The GGA potentials EV93 and AK13 as well as the gBJ 
potential are more anisotropic. The gBJ potentials leads 
to an improved agreement with EXX-OEP both in the 
bonding and interstitial regions, whereas the EV93 and 
AK13 potentials show too much variation in the inter¬ 
stitial region. We note that similar conclusions can be 
drawn for Si, BN, and MgO. 


However, it is also clear from Fig. [ 6 ] that the agree¬ 
ment between EXX-OEP and the best semilocal poten¬ 
tials (gBJ) is not perfect and that differences are still 
present. For a more detailed analysis we show in Fig. [7] 
one-dimensional potentials plots for Si. It becomes ev¬ 
ident from Fig. [7|^a) that AK13 (and to a lesser extent 
also EV93) leads to completely unphysical oscillations 
in the interstitial region of Si, while the EXX-OEP and 
gBJ potentials are rather flat and very similar to each 
other in this region. Actually, we have observed that in 
general the AK13 and EV93 potentials show such large 
oscillations in the wide interstitial regions present in such 
open structures, which is due to their enhancement fac¬ 
tors F x ( s ) in Eq. 0 whose magnitudes are much larger 
than for PBE and B 88 (see Fig. [l]). The direct effect of 
these much more positive values of the AK13 potential 
in the interstitial region is to shift up the unoccupied or¬ 
bitals (located mainly in the interstitial) relative to the 
occupied ones, thus explaining the positive (or less nega¬ 
tive than for most other potentials) AK13 values for the 
ME on the transition energies (Table [TT|). 

In Fig. gb) we can see that the height of the in¬ 
tershell peak at d 0.6 A is strongly underestimated 
and washed-out by PBE, whereas EV93 and AK13 tend 
to overestimate the peak height for Si. However, with 
increasing atomic number the ability of the AK13 and 
EV93 potentials to reproduce the height and position of 
the intershell peaks seems to improve. As shown in Fig.[ 8 j 
in the vicinity of the Cu atom both AK13 and EV93 
mimic the EXX-OEP quite accurately, while substantial 
differences are present at the O atom [see Fig.[ 8 |c)]. Sim¬ 
ilar observations hold for the intershells peaks in NiO. 
This is in agreement with Refs. [l7| and [28] which show 
that AK13 and EV93 reproduce very accurately the po¬ 
sition and height of the intershells peaks in transition- 
metal atoms. As already discussed in Sec. |IIB| and shown 
in Fig. [2| the height and position of the peaks with gBJ 
depend strongly on the three parameters 7 , c, and p. 
As shown in Fig. [7|[b), the intensity of the peak is too 
large with (q,c,p) = (1.4,1.1,0.50), but too weak with 
( 7 , c,p) = (0.6,1.0,0.60) (not shown). 

Concerning the BJ-based potential with the UC, gB- 
JUC, the results are very bad for the EXX total energy 
and energy position of core states as discussed above (see 
Tables [T] and [m|. This is a consequence of the very inac¬ 
curate gBJUC potential in the region close to the nuclei 


as shown in Figs.pfa) andpfc). As discussed in Sec. IIB 
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FIG. 8. (Color online) Exchange potentials v x in CU 2 O plotted from the Cu atom at site (1/2,1/2,0) (d = 0) in the direction 
of the O atom at site (3/4, 3/4, 3/4) (d — 3.54 A). Logarithmic scales on the x-axis were used for panels (a) and (c) which 
correspond to the vicinity of the Cu and O atoms, respectively. The potentials were shifted such that / cell v x (r)d 3 r — 0. 


the effect of replacing the kinetic-energy density t by D 
in the second term of Eq. (16) is very large in the core 


region of atoms with the consequence that the potential 
becomes too attractive. On the contrary, without the UC 
the gBJ potential resembles the EXX-OEP very closely 
in the core region [except at the position of the intershell 
peaks, see Figs. [SJa) and[ 8 ][c)]. 

Though, for CU 2 O it was mandatory to use the UC in 
order to obtain qualitative agreement with EXX-OEP 
for the transition energies and EFG. For the EFG in 
particular, this could seem puzzling that the gBJUC 
potential gives good results despite it looks quite inac¬ 
curate close to the Cu nucleus. However, as already 


mentioned in Sec. IV C the EFG is determined by the 
non-sphericity of the electron density p near the nucleus. 
More specifically, the EFG is determined essentially by 
the radial function plm with (L, M ) = ( 2 , 0 ) of the spher¬ 
ical harmonics expansion of the electron density inside 
the atomic sphere.® Figure [9](a) shows the (expected) 
very good agreement between the gBJUC, EXX-OEP, 
and HF methods for p 2 o (and also for p^o but not poo)- 
By looking at the corresponding radial function t > x?20 of 
the exchange potential [see Fig. [9][b)], we can observe a 
rather good agreement between gBJUC and EXX-OEP 
in the region beyond 0.2 A, which m ainly concerns the 
d-d component of the EFG (see Table |lV|) . We mention 
that the radial functions p 2 0 and v^o obtained with B 88 , 
EV93, and AK13 are qualitatively similar to PBE and 
gBJ without UC. For the p-p component, the agreement 
between gBJUC and EXX-OEP also comes from the va¬ 
lence region of the Cu atom and the interstitial, and, as 
shown in Fig. §b), these two potentials are relatively 
close to each other in this region. Actually, the correct 
description of the Cu-p states far away from the Cu nu¬ 
cleus affects the anisotropy of the Cu-p-states close to the 
Cu nucleus. These similarities observed in the EXX-OEP 
and gBJUC potentials can also explain the agreement for 


the transition energies. 

Turning to antiferromagnetic NiO, the difference rd — 
between the spin-up and spin-down exchange poten¬ 
tials is shown in Fig. [To) The angular e^-shape around 
the Ni atoms is the most pronounced with the EXX-OEP 
and gBJ [with ( 7 , c,p) = (0.4,1.3,0.65)] potentials, thus 
leading to the large band gaps between the t 2 ^ and e g 
states of the minority spin (see DOS in Fig. [5| in com¬ 
parison to the other potentials. However, it can also be 
observed that the magnitude of is the largest 

with EXX-OEP (i.e., compared to gBJ there are a cou¬ 
ple of additional isolines between the Ni and O atoms), 
which could explain the large exchange splitting between 
the spin-up and spin-down states observed in Fig. [5] for 
EXX-OEP. 


E. Analysis of the electron density 

Figure [TT] shows the electron density in Si obtained 
from various potentials minus the LDA density, which 
serves as reference. As discussed above for the case of C 
(Fig. §, which is very similar to Si and BN, the GGA, 
gBJ, and EXX-OEP potentials are more attractive (re¬ 
pulsive) than LDA in the bonding (interstitial) region 
of Si. Consequently, the electron density is increased 
(decreased) in the bonding (interstitial) region. From 
Fig. [Tl] it is rather clear that the gBJ and EXX-OEP po¬ 
tentials lead to very similar electron densities, while the 
EV93 and AK13 densities are, compared to EXX-OEP, 
too large (small) in the bonding (interstitial) regions. 

As shown in Fig. [ 12 ] for NiO, the effect of using a 
beyond-LDA potential is to increase the spin-up elec¬ 
tron density p^ on the Ni atom with a full spin-up 3d- 
shell (red regions around the Ni atom at the left up¬ 
per corner) and to decrease the spin-down density p± on 
the same Ni atom (which corresponds to p^ around the 
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EXX-OEP LDA 



FIG. 9. (Color online) Plots of the radial functions P 20 (a) 
and r » x ,20 (b) versus the distance from the Cu atom in CU 2 O. 
The functions are multiplied by r 2 . 


other Ni atom). This results in an increase of the mag¬ 
netic mom ent of the Ni atom as discussed above (see 
Table |lV|) . The other effect of using a beyond-LDA po¬ 
tential is to increase the ionicity (red regions around 
the O atoms). More quantitatively, compared to LDA 
the number of electrons inside the sphere of the Ni 
atom changes by +0.01 (PBE), —0.01 (EV93, AK13), 
-0.14 [gBJ(0.4,1.3,0.65)], and -0.11 (EXX-OEP), while 
for the O atom the changes are +0.10 (PBE), +0.16 
(EV93), +0.23 (AK13), +0.47 [gBJ(0.4,1.3,0.65)], and 
+0.35 (EXX-OEP), which shows that gBJ reproduces 
quite accurately the trends of EXX-OEP. From Fig.[l2j it 
is also rather clear that the gBJ(0.4,1.3, 0.65) and EXX- 
OEP electron densities are overall very similar (note in 
particular the asymmetry of the 3d density around the 
Ni atom with partially filled spin-up electrons). 

In Fig. [13] we show the density of the core electrons 
of the Cu atom in CU 2 O. Compared to the density p core 


FIG. 10. (Color online) Two-dimensional plots of the dif¬ 
ference between spin-up and spin-down exchange potentials 
(i>x — r’x) in a (001) plane of antiferromagnetic NiO. The con¬ 
tour lines start at —2 Ry (blue color) and end at 2 Ry (red 
color) with an interval of 0.235 Ry. The Ni atom with a full 
spin-up 3d-shell is at the left upper corner. 


obtained with EXX-OEP, the PBE core density is less 
contracted since it has smaller values for r < 0.14 A but 
is larger for r > 0.14 A. The reverse is true for the gBJUC 
potential since the positive values of pfori UC — 
are on average closer to the nucleus than the negative 
values, which is a consequence of the fact that close to 
the nuclei gBJUC is much more attractive than the EXX- 


OE P and all other potentials (as discussed in Secs. IIB 
and |IV D ). The AK13 and gBJ potentials lead to trends 
similar to PBE, however the discrepancies with respect 
to EXX-OEP are reduced. Note also that the gBJ poten¬ 
tial with the parameters ( 7 , c,p) = (0.6,1.0,0.60), which 
are more appropriate for the EXX total energy, leads to 
a slightly more accurate core density than with the val¬ 
ues (1.4,1.1, 0.50) that were determined for the transition 
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EXX-OEP-LDA PBE-LDA 



gBJ(0.6,1.0,0.60)—LDA gBJ(1.4,1.1,0.50)-LDA 



FIG. 11. (Color online) Electron density p obtained with 
different exchange potentials minus p LDA plotted in a (110) 
plane of Si. The contour lines start at —0.001 electron/bohr 3 
(blue color) and end at 0.001 electron/bohr 3 (red color) with 
an interval of 0.0002 electron/bohr 3 . 


energies. 


V. SUMMARY AND OUTLOOK 

In this work, we have compared several approximate 
semilocal exchange potentials to the exact EXX-OEP. 
The closeness between the semilocal and EXX-OEP po¬ 
tentials was quantified by considering the EXX total en¬ 
ergy and electronic band structure for various solids, as 
well as the EFG in CU 2 O and the magnetic moment in 
NiO. An attempt to parameterize a semilocal BJ-based 
potential has also been made and we have shown that 
by the introduction of a few parameters, it was possible 
to improve substantially the agreement with EXX-OEP 
compared to the GGA and original BJ potentials. How¬ 


EXX-OEP-LDA PBE-LDA 



mm 

EV93-LDA 

AK13-LDA 


# S' q 

gBJ(0.6,1.0, 0.60)—LDA 

gBJ(0.4,1.3,0.65)—LDA 




FIG. 12. (Color online) Spin-up electron density p^ ob¬ 
tained with different exchange potentials minus p^ DA plot¬ 
ted in a (001) plane of antiferromagnetic NiO. The contour 
lines start at —0.005 electron/bohr 3 (blue color) and end at 
0.005 electron/bohr 3 (red color) with an interval of 0.001 
electron/bohr 3 . The Ni atom with a full spin-up 3d-shell is 
at the left upper corner. 


ever, it became also obvious that there is no universal 
set of parameters that leads to satisfying results for all 
properties and solids at the same time. For instance, if 
a set of parameters is appropriate for the EXX total en¬ 
ergy (or electronic structure) of C, Si, BN, and MgO, 
then it will not work so well for antiferromagnetic NiO, 
and vice-versa. Another example was CU 2 O for which it 
is mandatory to use the UC to obtain qualitative agree¬ 
ment with EXX-OEP for the band gap and EFG, while 
the UC is very detrimental for the EXX total energy and 
energy position of the core states in all solids. 

From the results, it is clear but not surprising that 
although the BJ-based potentials lead to interesting re- 
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FIG. 13. (Color online) The difference between the density 
Pcore of the core electrons of Cu in CU 2 O calculated with 
semilocal functionals and p core from EXX-OEP (multiplied 
by 47rr 2 ). 


suits, the semilocal approximations show limitations and, 
furthermore, there is no systematic way to improve their 
accuracy. Beyond the semilocal level of theory, there is 
the group of exchange potentials which consist of the 
nonlocal Slater potential v ® [Eq. plus a term which 
is either nonlocal [like i n the K rieger-Li-Iafrate (KLI) 103 
or localized HF (LHF) 104 I 105 « potentials] or local with 
an event ual dependency on the energies of the occupied 
orbitals.®^® The computational cost of these poten¬ 
tials is rather high since the Slater potential and the 
nonlocal terms require the calculation of HF-type inte¬ 
grals (see Ref. [83] for a summary of the expression of 
these potentials). These nonlocal potentials avoid some 
of the technical difficulties of EXX-OEP as their con¬ 
struction involves only the occupied orbitals. There are 
numerous studies on the Slater-based potentials and we 
just mention Ref. 11071 where it was shown that the KLI- 
and LHF-generated orbitals lead to EXX total energy of 
atoms which are much lower than with BJ orbitals. How¬ 
ever, Engel (Ref. [60]) noted that the KLI approximation is 
not able to open the band gap in antiferromagnetic FeO, 
while a band gap of 1.66 eV is obtained with EXX-OEP 
(augmented by LDA for correlation). 


On the other hand, as already mentioned in Sec. |IIB[ 
the semilocal BR potential— [Eq. (10)] seems to repro¬ 
duce (at least visually) quite accu rately the features of 
the Slater potential in atomsf 22 ^ however, the agree¬ 
ment is not perfect (see Ref. [33)) and the comparison of 
the EXX total energies evaluated with BJ(Slater) and 
BJ(BR) orbitals shows non-negligible differences in some 
casesP 2 Avoiding the calculation of the Slater potential 
by using the BR potential instead would certainly be ad¬ 
vantageous, but more comparison studies between the 
BR and Slater potentials are needed. 


A possible way of improving the reliability of a semilo¬ 
cal potential (e.g., gBJ) to reproduce EXX-OEP, could 
be to use a similar parameterization as the one used for 
the constant c in the mBJ potential [Eq. (15)]: 



where a and /3 are parameters and V^ e u is the unit cell 
volume. It has been shown that with the optimized val¬ 
ues a = —0.012 and /? = 1.023 bohr 1 / 2 , mBJ reproduces 
with rather great acc uracy the experimental band gap of 
many solids IMI 80 l 108 r 10 9 Actually, the use of an integral 
expression like Eq. © is a way to introduce nonlocal¬ 
ity (similar as with the Slater potential), but in a cheap 
way since there is no summation over orbitals like in the 
Slater potential. However, the drawback of using the av¬ 
erage of |Vp| /p in the unit cell is that this quantity is 
infinite for systems with an infinite vacuum (e.g., isolated 
molecule or surface). An alternative to Eq. © which 
can be applied to any kind of systems might be helpful in 
improving the universal character of a potential like gBJ. 
For instance, as suggested by Marques et al in Ref. 11101 
a possibility would be to make c r-dependent and the 
integrand in Eq. © localized around r by multiplying 
|Vp|/p by a function of |r — r'| which goes to zero at 
|r — r'| 00. 

In order to adjust the parameters of an approximate 
functional for each system, an approach as suggested in 
Ref. 11111 might be helpful. The free parameters of the 
potential are adapted at each iteration such that the 
EXX total energy becomes minimal. However, such a 
procedure is rather expensive since the equations to de¬ 
termine the parameters involve HF-like matrix elements 
between occupied and unoccupied orbitals. Neverthe¬ 
less, this would be a way to adjust the parameters for 
each solid and therefore improve the universality of the 
potential. 

More over , we mention the work of Star over ov and co- 
worker ^ 3 * 107 ! who proposed an expression for a multi¬ 
plicative exchange potential (making no use of unoc¬ 
cupied orbitals) which leads to results that are quasi¬ 
identical to the EXX-OEP results. However, this ap¬ 
proach requires the HF orbitals which reduces its use for 
solids and large scale applications. 

Finally, we note the very few studies reporting OEP 
calculations including correlation like the random-phase 
approximation (RPA) in addition to EXX (see Refs. [47j 
I4H1 [5811591 and l6H for results on solids). The RPA-OEP 
potentials could certainly also serve as reference for the 
modelling of realistic multiplicative exchange-correlation 
potentials u xc including correlation. 
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